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SECTION 1.0 
INTRODUCTION 


This report describes the program called CTRANS which 
has been designed to perform radiative transfer computations 
in an atmosphere with horizontal inhomogeneities (clouds) . 
Since the atmosphere-ground system was to be richly detailed, 
the Monte Carlo method was employed. This means that results 
are obtained through direct modeling of the physical process 
of radiative transport. The effects of atmospheric or ground 
albedo pattern detail are essentially built up from their 
inpact upon the transport of individual photons. CTRANS 
actually tracks the photons backwards through the atmosphere , 
initiating them at a receiver and following them backwards 
along their path to the sun as shown schematically in 
Figure 1.0.1. Backwards tracking has several advantages: 
the pattern of incident photons generated through backwards 
tracking automatically ref lects the importance to the receiver 
of each region of the sky. Further, through backwards 
tracking, the impact of the finite field of view of the 
receiver and variations in its response over the field of 
view can be directly simulated. The backwards tracking 
method, additionally, is well suited for segregating results 
according to the character of the final scattering suffered 
prior to entering the receiver. 

Section 2.0 describes the backwards tracking Monte 
Carlo method. In Section 3.0, the results of a range of 
test applications are reported which serve to illuminate 
the capabilities and limitations of the program. Section 

4.0 contains brief descriptions of each routine comprising 
CTRANS. Since CTRANS is a highly generalized program with 
many operational modes and a variety of data structures, 
Section 5.0 outlines the input parameters. Finally, Section 

6.0 contains a FORTRAN listing of CTRANS. 
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TRACK OF HYPOTHETICAL 
PHOTON WHICH IS SAMPLED 



GROUND PLANE 


FIGURE 1.0.1. SCHEMATIC REPRESENTATION OF THE PROBLEM GEOMETRY 
AND THE BACKWARDS MONTE CARLO METHOD 
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SECTION 2.0 
MONTE CARLO METHOD 


CTRANS employs the Monte Carlo method to compute flux 
or intensity upon a receiver situated in or above a plane 
parallel atmosphere. Photons are tracked backwards through 
the atmosphere from the receiver along their path to the sun. 

The particular advantage offered by the Monte Carlo 
method lies in the degree of complexity of the earth- atmospheric 
system which can be modeled. 

The present code can handle a plane parallel atmosphere 
with up to 100 homogeneous layers. Each layer can be composed 
of up to five distinct scattering (or absorbing) species: 
molecular scatterers, monodisperse or polydisperse aerosols 
or hazes, ozone, etc. Further, there may be up to 10 finite 
rectangular solid clouds which are superimposed upon the 
layered atmosphere. Each cloud may be composed of any of the 
five allowed active atmospheric components with an independently 
specified density for each cloud. Thus, the atmosphere may 
be inhomogeneous both vertically and horizontally. 

The treatment of ground reflection is also generalized. 
Reflectance at the ground may be either Lambertian or Fresnel 
in character. Fresnel regions are further characterized by 
a surface roughness specified by a probabilistic slope distri- 
bution (so chosen as to provide an approximate model of a 
sea surface with waves). The ground plane need not be uniform; 
rough Fresnel and Lambertian regions may be intermixed in a 
patchwork fashion as specified by a reflectance type map. 

Any Lambertian regions may be further broken down into patches 
having different Lambertian albedos. Results are computed 
for up to 50 different overlying maps of the Lambertian albedo, 
simultaneously. Aside from the gain in computational efficiency, 


3 



this is important because it provides a means for computing 
the relative effects of changing the ground albedo pattern 
with a precision which far surpasses the absolute precision 
of the Monte Carlo results for a single map. 


2.1 PHOTON TRANSPORT AND SAMPLING 

CTRANS tracks photons backwards initiating them at the 
receiver and following them through successive scatterings on 
their path through the atmosphere. Backwards tracking has been 
adopted for a number of reasons. It permits the efficient and 
direct simulation of receivers with finite fields of view and/ 
or finite extent. A forward tracking scheme would require 
the distribution of initial photons over a wide area when 
horizontal inhomogeneities are present and this would introduce 
large variances in the results. By initiating photons at 
the receiver, this source of variance is largely avoided. 


To understand backward tracking, consider first a 
single scattering event. If 1^ is the Stokes vector of 
incident unscattered light, S is the scattering phase matrix 
and R is a rotation matrix which serves to rotate the plane 
of reference of the Stokes vector into the scattering plane, 
the Stokes vector after scattering I_^ is given by 



R 


SRI 
o 


where R rotates the final Stokes vector into the reference 

T 

plane of the receiver. After two scatterings, the final 
Stokes vector will be 



Cast in this way, it is apparent that what is important is 
the product matrix which is built from the scattering matrices 
appropriate to the string of scatterings suffered by the photon 
This product can be built up either from the left or from the 
right with the same result; i.e., we may equally well follow 
the photon’s history forward or backwards if we accumulate 
scatterings properly into the product [which we call the 
cumulative scattering matrix] . 

There is a further advantage in tracking backwards: 
contributions from several solar directions may be accumulated 
from a single track. For example, if the photon is tracked 
backwards from the receiver, only the last rotation and 
scattering contain any reference to solar direction. Replacing 
one sun by another means changing only these two matrices. 

For each solar direction, at most, two 4x4 matrix multipli- 
cations must be performed. 

Photons are initialized at the reciever; which may be 
general in the following ways: it may have finite or infinites 

simal (conical] field of view; it may be oriented pointing 
toward any zenith angle with any azimuthal angle; it may be 
at any altitude and have any specified x-y coordinate. On 
the other hand, for flux calculations the receiver aperture 
may be finite. This is accomplished by letting the receiver 
be an entire cloud face (any face of any cloud, pointing into 
or out of the cloud) . 

First, initial photon coordinates are created. If the 
receiver is a cloud face, initial photon coordinates are chosen 
from a uniform distribution over the cloud face. The photon 
direction is then determined using: 1) uniform distribution 
over solid angle for a cloud-face receiver, or 2) cosine distri 
bution over the conical field of view of an inf initessimal 
receiver, or 3) in a specified direction if the field of view 
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is infinitess imal . For each type of distribution, the appro- 
priate weight is computed (to be multiplied into the sampled 
values upon sampling). 

The distance to be traversed to a scattering is computed 
on the basis of optical distance to be traversed: if p is a 
uniformly distributed random number, the optical distance to 
be traversed is given by 

x = - tap (2.1.1) 

Through reference to a table containing the optical thickness 
of each layer of the ambient atmospher and taking into account 
the added optical density in the interior of clouds, the optical 
distance traversed and physical distance traversed are incre- 
mented until an optical distance r has been reached. If, before 
this has occured, the photon exits the atmosphere, the photon 
is terminated. If the photon strikes the ground, the scattering 
will be a ground scattering and a ground scattering indicator 
is set (ISCAT=1) . 

Having located the position of a scattering, the scatterer 
type is chosen on the basis of the relative contribution of 
each scattering species at that point to the local optical 
density. Here, a cumulative distribution function, which has 
been tabulated for each layer of the ambient atmosphere, is 
used in conjunction with any additional optical density due 
to the presence of a cloud. Photons are not explicitly absorbed. 
Instead, a weight is accumulated which at any point along the 
path reflects the probability that the photon has not been 
absorbed. 

After having determined the scatterer type, the probability 
that the photon could have come from each sun is sampled. If 
Sj is the cumulative scattering matrix at the point of sampling, 
the sampled Stokes vector will be 
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I 


( 2 . 1 . 2 ) 




R. S. 


_2_ _l 



- 1 ■; 
e J 


w 


Where R- is a rotation matrix appropriate to a scattering into 

tji ^ 

the j solar direction; Sj is a scattering matrix (a Mueller 
matrix if the scatterer is an atmospheric constituent) ; I is 
the incident Stokes vector; t^ is the total optical distance 
along. the path connecting the scattering point and sun^ ; w 
is a weight reflecting the probability that the photon has 
not been absorbed (including ground absorption for Lambertian 
regions encountered). There will be a different value for 
w for each Lambertian albedo map. Note that the above form 
assumes that 1^ represents intially unpolarized light. If 
the incident light from the sun were polarized, an additional 
rotation would be necessary. Rj has the form 



I 1 ° 

/ 0 cos 2X 

I 0 -sin2X 
\ 0 0 


0 

sin2X 

cos2X 

0 



(2.1.3) 


where X is the rotation angle between the normal to the pre- 
vious scattering plane (or the reference direction in the 
case of first scatterings) and the normal to the scattering 
plane for a scattering into the solar direction. has 

the general form 




(2.1.4) 


where S„ m 
£ ,m 

scattering , 


depends upon scattering angle. 

s n - S 22 - j a ♦ c °= 2 V 

S 12 = f ( cos2 6-1) 


For Rayleigh 
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S 33 S 44 ” 7 cos 


S 34 - 0 


where 3- is the scattering angle for a scattering into the 
th J 

j solar direction. For Lambertian scattering at the ground, 

th 

S, , = 4 cos 9., where 8. is the zenith angle of the j sun 
11 3 3 

and all other matrix elements are zero. Lambertian reflection 
thus is perfectly depolarizing. 


Mie scattering and Fresnel scattering both have matrices 
with the form of Eq. (2.1.4) . Mie elements are stored in 
tabular form and values for particular scattering angles are 
obtained by interpolation. The Fresnel matrix elements are 
given in the discussion of ground reflectance. 


The next step after sampling is to accumulate a "real" 
scattering. Knowing the scatterer type, an appropriate scatter- 
ing angle is found. If the scattering angle were chosen from 
a uniform distribution, and if the scattering matrix appro- 
priate to the medium exhibited sharp peaks, then most photons 
would rapidly lose weight and large numbers of photon tracks 
would be necessary in order to achieve a reasonably small 
variance. To resolve this difficulty, we use the concept of 
importance sampling. We choose the scattering direction 
from a biased distribution designed to maximize the resultant 
contribution. The bias thus introduced is then removed by 
applying a weight. More concretely, we wish to fold in a 
scattering matrix 
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0 


0 


1 

f(3) 




S 11 S 12 

S 12 S 22 
0 0 

0 0 


0 0 
S 33 S 34 

" S 34 S 44 


) 


f(B) dw 
4u 


(2.1.5) 


Instead of choosing the scattering angle from a probability 
distribution corresponding to the probability density 
(uniform), we use the distribution corresponding to the 
probability density — The weight of the photon follow- 
ing the scattering can thus be maximized by choosing f(3)=S 11 
times a constant chosen so that the first matrix element is 
always unity. A real scattering is accumulated by multiplying 
into the cumulative scattering matrix first a rotation, and 
then the scattering matrix. 

The process is repeated until either a designated max- 
imum number of scatterings have been suffered, the photon 
leaves the atmosphere or its weight drops below the assigned 
threshold (due, for example, to atmospheric absorption). 
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2.2 


ENCOUNTERS WITH CLOUDS 


The presence of finite clouds complicates the evaluation 
of the fate of the path of the photon in that there is no simple 
way to determine if a given straight line in space enters a 
specified rectangular box. It was therefore necessary to 
undertake an analysis to find optimum procedures for evaluating 
the path as quickly as possible. 

This analysis is applied in two areas of the program: 

1. In evaluating the path of the ’’backward” photon 
to the "previous" scattering. 

2. In computing the probability that a "forward" 
hypothetical photon would reach a given scattering 
point from the Sun without scattering. 

The problem has been treated in two parts. First, we 
have found a quick and efficient method for eliminating clouds 
which could not possibly be entered by the photon path, and 
secondly, we have performed a detailed analysis to evaluate 
the uneliminated clouds. 

There is a simple method for determining whether it is 

possible for a straight line to enter a given rectangular cloud. 

The approach is to circumscribe the cloud with a sphere, the 

center of which coincides with the center of mass of the cloud. 

2 2 21/2 

The radius, r, of the sphere can be taken to be (a +b +c ) ' 
where the lengths of the sides of the cloud are (2a, 2b, 2c). 

In this case the eight corners of the box lie on the sphere 
as illustrated in Figure 2.2.1. 
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Since the sphere completely contains the cloud then 
if a line cannot enter the sphere, it cannot enter the cloud. 
Thus, the first check should be to find if the path can enter 
any cloud circumscribing spheres. We shall present the 
analysis for a single cloud with the understanding that this 
must be repeated for all clouds in the program. 

Suppose we have a line with direction cosines (c ,c ,c ) 

x y z 

leaving a point A(x^,y^,z^). The equation of the line is 


x-x, y-y-i z-z. 


( 2 . 2 . 1 ) 


The distance, R, of the point 0 (X,Y,Z) from A is 


R = { (X-x . ) L + (Y-y.r + (Z-ZtT) 


1/2 


( 2 . 2 . 2 ) 


and the point of closest approach of the line to 0 is such 
that AP is perpendicular to OP, as illustrated in Figure 2.2.2. 
Now the cosine of the angle OAP is given by the dot product 
of the direction cosines of AP and AO; i.e.. 


cos OAP 


(X-X P C x + (Y-yPSr + < Z - z lK 
R R R 


(2.2.3) 


But the distance AP is given by R cos OAP; i.e.. 


AP = (X-x 1 )c x + (Y-y 1 )c y + (Z-z 1 )c z 


(2.2.4) 
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and by Pythagorus's theorem, 


2 2 2 
OP = R - AP 

= R 2 - {(X-x 1 )c x + (Y-y 1 )c y + (Z-z 1 )c z ) 2 


(2.2.5) 


2 2 

Now, if OP is greater than r , the squared radius of the 
circumscribing sphere, then the path cannot enter the sphere. 
Also, if the quantity in curly parentheses is negative, then 
again entry is impossible. 

2 

Note that OP can be computed in the general coordinate 
system and that no square roots or complicated functions are 
necessary for its evaluation. 

If we are tracking the path of the backward photon then 
another check can be made. Although the point of closest 
approach may be inside the circumscribing sphere, it may be 
that the length, L, of the path is insufficient to reach the 
sphere. Suppose, in Figure 2.2,2 that the line reaches the 
sphere at Q. Then OQ has length r, the radius of the sphere. 
Now from (2.2.4) and (2.2.5) we are able to evaluate both 
OP and AP and by Pythagorus ' s theorem we have that 


QP 


2 




( 2 . 2 . 6 ) 


Thus we are able to evaluate the distance to the sphere as 


AQ = AP - QP 


(2.2.7) 
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which may be evaluated by (2.2.4) and (2.2.6). Testing this 
value against L gives us another possible means of cloud 
elimination . 


We have attempted to use the circumscribing sphere 
concept further, but we finally concluded that there are 
some dangers in its further application. For example, the 
sequence of distances along the path to a cloud need not be 
the same as the sequence of distances to the circumscribing 
spheres or the point of closest approach. An illustration 
of such a case is given in Figure 2.2.3. The reason for the 
difficulty is that the circumscribing sphere gives no infor- 
mation as to cloud orientation and can therefore mask orien- 
tation differences which may be important. For this reason 
we proceed in the program at this point to a completely 
general analysis. 

The first step in treating an uneliminated cloud is 
to convert to the coordinate system of the cloud. Thus the 
location of the scattering point (x^y^z^) is transformed 
t° ( x i c >yi c >z lc ) using 

X lc = ( X i' x ) COS0 C + ^i“ Y ^ sin0 c 


y lc = “( x i _ X) sin© c + (y^Y) cos9 c (2.2.8 



where 0 c is the orientation angle of the cloud and (X,Y,Z) 
are the coordinates of the center of mass of the cloud. The 
direction cosines of the line are transformed to the cloud 
system using 
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FIGURE 2.2.3. EXAMPLE OF SEQUENCING FAILURE OF 
CIRCUMSCRIBING SPHERES ALGORITHM 
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C cos0 + C sin0 
x c y c 


cx 


C =-C sinO - C cos0 
cy x c y c 


C = C 
cz z 


(2.2.9) 


We shall present the analysis for determining if the 
photon enters a side perpendicular to the cloud x axis. In 
the computer program we index terms relating to the cloud 
x, y and z axes 1, 2 and 3, respectively. The analysis for 
the sides perpendicular to the y and the z axes can then be 
performed by merely rotating the ordered index set (1,2,3) 
to (2,3,1) and (3,1,2), respectively. 

The following tests will be performed in sequence. 

1. Is l x 2 C l £ a, the half length of the edges (2.2.10) 

parallel to the cloud x axis. 

If this condition is satisfied, then the photon 
cannot possibly enter a face perpendicular to the 
cloud x axis. 

2. If x^ c < -a, then the path can enter the face 

x =-a. In this case, set q=-a. 
c n 

If x-^ c > a, then the path can enter the face, x c =a. 

In this case, set q=a. 

3. Compute (x. -q) *C . If this is positive, then the 

XL C a 

photon is traveling away from the face, x c =q, and 
will not enter the cloud. 
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4. Compute the (y c , z ) values of the point of inter- 

section of the line with the plane, x c =q. These 
are given by 


y c i = y ic + c 


q-x 


lc 


cy 


cx 


q-x 


lc 


z T = z, + C 
cl lc cz 


cx 


( 2 . 2 . 11 ) 


Check to see if both the conditions (2.2.12) apply. 


y cl I < b 


z cl' < C 


( 2 . 2 . 12 ) 


If these are both true, then the photon enters the 
cloud at the point (q, y cI , z c j) • Otherwise, entry 
does not occur at the x=q plane. 

When it is determined that a photon has entered a 
given cloud then the tracking mechanism from that point on 
will be performed in the cloud coordinate system. 

If it is determined that a later scattering point 
(x lc ,y lc , z lc ) lies outside the cloud then we determine the 
exit point using steps (4) and (5) only. 
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For evaluating the track of the "forward" hypothetical 

photon from the Sun to a scattering center we will determine 

if the direction cosine C of the line from the scattering 

cxs ° 

point to the Sun is positive or negative. (The solar direc- 
tion cosines (C , C , C ) will be computed for each 
cloud prior to the simulation.) If C cxs is positive, we set 
q=a, whereas if it is negative, we shall set q=-a. At this 
point steps (4) and (4) are employed to find the point of 
exit from the cloud. 
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2.3 


GROUND SCATTERING 


CTRANS has the capability to model reflections from 
a spatially inhomogeneous ground. In the present version 
of the code, the physical mechanism of reflection may be 
modeled as either Lambertian reflection or as Fresnel 
reflection from a rough surface and these modes may be 
intermixed in a patchwork fashion. 

Reflectance type is basically controlled by a reflec- 
tance type map which dictates the reflectance character from 
point to point across the infinite ground plane. When a 
photon strikes the ground, its impact coordinates are com- 
puted and then used to decide whether the point of impact 
was within a Lambertian or a rough Fresnel region. 

Since Lambertian reflection is not physically depen- 
dent upon the direction of incidence or reflection but only 
upon a scalar Lambertian albedo, CTRANS can model simul- 
taneously the effect of several values (up to 50) of the 
albedo for each of the Lambertian regions. This is accom- 
plished in the following way: A vector of scalar weights 

is maintained; one for each albedo map. When a Lambertian 
scattering occurs, the albedo at the point of impact is 
determined for each of the albedo maps and this is multiplied 
into the corresponding albedo weight to obtain a new weight 
for subsequent use. Upon sampling, each current albedo 
weight is folded into the sampled Stokes vector to obtain 
samples corresponding to each of the allowed albedo maps. 

If the impact point lies within a Fresnel region, however, 
the vector of albedo weights is left unchanged (since it 
represents the cumulative effects of only the Lambertian 
regions) . 
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Symbolically, the entire reflection process is as 
follows: given an old cumulative scattering matrix S 0 , the 

new one, S^, is computed as: 


S,, = S R s 

N o 


where R is a rotation matrix and s depends upon the reflec- 
tance type at the impact point: 


Lambertian s 


Fresnel s = 



where 



represent Fresnel matrix elements: 


11 

(R n 2 * 

R 1 2 )/2 

12 

CRu 2 - 

R x 2 )/2 

33 

R 11 R 1 
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and 



n cos x t " cos x t 
n cos x jL + cos X t . 


cos X-: - n cos x t 

r i = -1 1 

cos \ i + n cos x t 

n = relative index of refraction of the surface 

X^ = angle of incidence 

X t = angle of refraction, computed from x- through 

Snell's law at the interface: n sin x t ~ sin • 



X,Y = coordinates of impact point 

J=0 = region is Lambertian and the values of 

a^(i=l,50) are determined from the Lambertian 
albedo maps according to the coordinates CX,Y) 

J= 1 = region is Fresnel a^ = l(i=l,50). 
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The two kinds of ground maps (Lambertian albedo maps 
and the reflectance type map) are specified in CTRANS in 
two different ways. 

The Lambertian albedo maps are specified in a fixed 
form by the subroutine ALBEDO. To change albedo map struc- 
tures one version of ALBEDO is interchanged with another. 

The possible variety of maps is thus limited only by the 
imagination. ALBEDO, in addition to its function of provid- 
ing specific values of the scalar weights at a specified 
point has the function of writing out text to describe the 
various maps. 

The reflectance type map is specified by options read 
in by subroutine INPUT and is therefore less flexible than 
the albedo maps. As options, the following choices are 
presently available: 

1) Uniform plane. 

2) Adjoining half planes of different reflectance 
type bounded by a line parallel to the y-axis 
whose position is specifiable by its x- coordinate . 

3) Non-overlapping rectangles set upon a uniform 
background (up to 20 such rectangles are allowed) . 

In all of the above, the background is specifiable as either 
Lambertian or rough Fresnel with inhomogeneities having either 
rough Fresnel or Lambertian character, respectively. The 
rectangles are specified by the (X,Y) coordinates of their 
centers and their X-Y edge lengths. 
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Text describing the chosen pattern is written out by 
INPUT. A cautionary note is in order here. Since the two 
types of maps are specified by unrelated means, care must be 
taken in the interpretation of results: The Albedo maps 
only apply to regions specified to have Lambertian character. 
Thus, using 50 varigated albedo maps and specifying the 
reflectance type map as a uniform rough Fresnel plane would 
produce identical results for all 50 albedo maps. 

For rough Fresnel regions, the roughness is characterized 

as a distribution of randomly oriented sloping plane facets. 

The probability density for the distribution of slopes was 

chosen to have two-dimensional Gaussian form. p(Z , Z )x 

x y 

6Z 6Z is the fraction of a small horizontal unit area of 
x y 

surface for which the x,y components of the slope are within 

the limits Z + i SZ and Z + i 6Z . We specify p(Z , Z') as: 
x z x . y z y xy 

P(z x , z y ) - (TO 2 )' 1 exp t-(Z x 2 * z y 2 ) / o 2 ] 

2 

The variance, a is specified either by three parameters WO, 

W1 , W; or as a set value specified on input. When specified 

2 

through parameters, a is given as 


o 2 = WO + W1 x W 


W represents wind speed (m/sec) in a model for the surface of 
the sea 1 -* and WO, W1 have either their default values (WO = 
.0015, W1 = 2.54 x 10 ) or values specified by card input. 

In the present code, all rough Fresnel regions share the 
same parametric values. This restriction could be relatively 
easily relaxed if desired. 
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The above ^form for p(Z , Z ) represents that applicable 

x y 

for a photon incident from the vertical. For photons incident 
from some other direction, the distribution must be modified. 

/v 

Let the photon's direction be given by unit vector k, the 
direction of the normal to a facet be n, and the vertical 

A 

direction be Z. The probability of encounter will then be 
given by: 


k*n 

-s-s — P(n) . 
(n-Z)Ck-Z) 


Having determined a slope at the photon impact point, 
the angle of incidence and, thence, the Fresnel matrix 

t 

elements can be determined. If k and k are the directions 
of propagation of the incident and reflected photons, their 
relation to n is 


A | A 



When sampling the contributions of unscattered photons from 
the Sun at a Fresnel scattering point, the direction of 
reflection is given by the solar direction and one merely 
evaluates the probability of having encountered a slope 
having the proper orientation. Note that in the present 
model the location of the Sun is given by a Dirac delta 
function (incoming plane waves) rather than being spread 
over a small angular region. The implied integration over 
the solar direction (or, equivalently, the change of vari- 
ables within the delta function specifying the solar 
direction) introduces an additional factor of 

A A 7 A _ 1 

.25 (n-Z)"- 5 (k-n) 
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Our present formulation does not properly account for 
shadowing (waves whose surfaces are inaccessible due to their 
orientation) or for hiding (slopes inaccessible because of 
intervening higher waves). Thus, we expect some systematic 
bias, especially when either the sun or the receiver is low 
in the sky. 
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SECTION 3.0 

TESTING: CAPABILITIES AND LIMITATIONS 

OF CTRANS 


CTRANS has been tested and validated through compari- 
son with the results of other calculations under a variety 
of conditions. Here, the results of representative cases 
will be presented and discussed. 

The basic functional structure of the Monte Carlo code 
(backwards-tracking algorithm, Lambertian ground reflection, 
statistical normalizations, etc.) may be tested via calcula- 
tions of radiative transfer through a pure Rayleigh atmos- 
phere and comparison with the well established results of 
Coulson, Dave, and Sekera. } Next, the atmosphere is made 
more complicated by adding layered structure and additional 
active atmospheric species. This provides a basic test of 
the algorithms for scatterer selection, the treatment of 
absorption, and the treatment of polydispersions. There may 
be special problems associated with optically thick atmos- 
pheres composed of strongly anisotropic scatterers so this 
case is of interest. 

CTRANS' capacity to handle horizontal inhomogeneities 
is basically evaluated by means of a calculation with a 
single cloud in the atmosphere and comparison with the 
results of McKee and Cox. } In this test, some difficul- 
ties were detected which are associated with the calculation 
of flux (in optically thick clouds) by means of a backwards 
tracking code. These difficulties are partially resolved by 
using a specially adopted forwards - tracking code which will 
be briefly described. 
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CTRANS' capacities to treat horizontal ground reflec- 
tion inhomogeneities are presented through test results. 
Absolute validation is difficult here in the absence of 
dependable analytical results, but having validated the 
microscopic reflection mechanism (through the inspection of 
results for uniform ground planes) only problems associated 
with the statistical density of sampling might be present. 
These statistical problems may be detected through the 
comparison of results from statistically independent runs. 


3.1 HORIZONTALLY HOMOGENEOUS ATMOSPHERES 

The first and most basic test case is that for a 
horizontally (and vertically) homogeneous Rayleigh atmosphere. 
Here CTRANS results may be compared with those of Coulson, 

Dave and Sekera J who have tabulated values for the Stokes 
parameters and polarization for Rayleigh atmospheres with 
optical thicknesses ranging up to 1.0 and for uniform 
Lambertian ground planes \\fith albedos from 0.0 to 0.8. 

We have run tests over a range of optical thicknesses. 

Of these, the most stringent test is provided by the largest 
optical thickness (t = 1). Since, in this case, single 
scattering, multiple scattering and ground reflection all 
are relatively important. Typical results for I, Q, U, 
and the polarization are shown in Figures 3.1.1 through 3.1.4. 
These results are for a receiver on the ground looking up with 
a zenith cosine angle of 0.4, 60° from the solar- vertical plane. 
As can be seen, agreement is excellent both for zero ground 
albedo and high ground albedo (0.8). Thus reflections at 
the ground are handled properly as is the basic transport 
through the atmosphere. 
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FIGURE 3 . 1 . 1 . 

I VS. SOLAR ZENITH COSINE 
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(SOLAR ZENITH COSINE) 


FIGURE 3.1.2. 

Q VS. SOLAR ZENITH COSINE 
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(SOLAR ZENITH COSINE) 




(SOLAR ZENITH COSINE) 




FIGURE 3.1.4. 

P VS. SOLAR ZENITH COSINE 
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(SOLAR ZENITH COSINE) 



In the next test, we add vertical structure and a 
polydisperse aerosol in a fairly realistic model of a clear 
atmosphere. Since analytical results are unavailable for 
such an atmosphere, we compare CTRANS' computed values with 
two sets of results kindly provided by Dr. R.S. Fraser 
(computed by programs named RADTRAN and VPD) . Results are 
presented in Tables 3.1.1 and 3.1.2 for the intensity and 
polarization at the bottom and the top of the atmosphere. 

As may be seen, CTRANS is in agreement with both routines 
to within the bounds expected from CTRANS' computed standard 
errors . 


CTRANS may encounter difficulty in handling optically 
thick atmospheres composed of scatterers having an extremely 
peaked phase function such as is provided by a Deirmandjian 
Cl polydispersion. The manifestation of the problem is a 
rather large variance especially evident when computing the 
radiance for a receiver looking down into the atmosphere. 

Tables 3.1.3 and 3.1.4 provide an illustration of the effect. 
Here we present a test case for a receiver pointing towards 
the nadir from the top of the atmosphere. The atmosphere 
is composed of pure Mie scatterers (Cl Deirmendjian haze: 
a = 6.0, B = 1.5, y = 1.0, index of refraction = 1.5). The 
wavelength is 0.825 microns. In the tables, I c refers to 
CTRANS' results; 1^. are results computed by a doubling 
method (provided by Mr. Lee Kyle); and SE c denotes the stan- 
dard error inherent in I . For this case, the ground albedo 
was zero and polarization was neglected. As may be seen, the 
agreement between I and 1^. for the x = 0.1 case is acceptable. 
However, when the optical thickness is increased to x = 10.0, 
the variance rises markedly (though no bias is detected) . 

This variance is extremely difficult to reduce; particularly 
for the x = 10.0 case. Here 2000 photons were tracked and 
up to 50 scatterings allowed per track, resulting in 39,422 
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Table 3.1.1. Realistic Clear Atmosphere Test: Intensity 
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Table 3.1.2. Realistic Clear Atmosphere Test: Polarization 
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Table 3.1.3. Uniform Aerosol Test t=0.10 
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scattering events (19.71 scatterings per photon track on 
the average). Tracking 2000 photons in this thick atmos- 
phere required 5.10 minutes of CPU time. The variance might 
be reduced by increasing the number of photon tracks, but in 
general this would be very expensive. There are two effects 
operating in the x = 10. Case: 1) Since the atmosphere is 

thick, photons must scatter many times before wandering out 
of the atmosphere. 2) The fact that the scattering phase 
function is extremely peaked in the forward direction means 
that most photons will not contribute importantly to the 
sample; i.e., only a few photons which, by chance, are travel 
ing nearly toward the sampling direction provide large con- 
tributions to the sample while all others contribute small 
amounts. The backwards tracking code is thus statistically 
inefficient in this case. This inefficiency may be unimpor- 
tant in optically thin atmospheres (or in absorbing atmos- 
pheres) since it may be economically overcome by increasing 
the number of photon tracks. 

In an attempt to improve the efficiency of the code, 
a number of alternate importance sampling schemes were imple- 
mented. In all of these, it was found that any substantial 
reduction in variance was accompanied by the introduction of 
bias. These alternate schemes have, thus, been abandoned. 

3.2 HORIZONTAL ATMOSPHERIC INHOMOGENEITIES (CLOUDS) 

In order to validate CTRANS' capacity to handle finite 

3) 

clouds, we may compare with the results of McKee and Cox. ' 
Their computation utilized a Monte Carlo method to calculate 
the flux emerging from the faces of a finite cloud of cubic 
shape. Their model ignored scattering outside the cloud and 
neglected polarization effects. Accordingly, we used no 
ambient atmosphere and restricted the ground to zero albedo. 
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McKee and Cox obtained results for clouds composed of an 
aerosol having a Deirmendjian haze Cl^ size distribution 
and optical thickness per kilometer equal to 4.9 and 73.5. 

The wavelength is 0.45 ym. 

The scattering phase function for this case contains 
a very narrow, strong forward peak. Our previous considera- 
tions would lead us to expect that these clouds would provide 
a very severe test of CTRANS. Optically thick clouds entail 
many scatterings per photon track and a sharply peaked phase 
function requires many independent photon tracks to achieve 
satisfactory variance. On this account, we have restricted 
the test to a cloud with an optical thickness of 4.9 per 
kilometer . 

The basic geometry of the test is illustrated in 

Figure 3.2.1. Light is incident upon the cloud from zenith 

angles of 0°, 30°, and 60°. The incident flux is normalized 
2 

to 1 watt/km normal to the beam. 

In order to compare with McKee and Cox, CTRANS (using 
backward tracking) initiates photons from a chosen face of 
the cloud, distributing them uniformly over the face and 
distributing their initial flight directions over 2 tt stera- 
dians. In so doing, CTRANS might be expected to incur addi- 
tional variance since some regions of the cloud will be 
brighter than others and some exit directions will be more 
important. This, then, will represent an additional source 
of difficulty for the backwards tracking method. For this 
test, separate computations are performed for each face. 

This means that flux will not be automatically conserved and 
the degree to which flux is conserved will be a test of the 
program . 
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When photons are initiated on a face, they are segre- 
gated according to whether the (forward) photon would be 
emerging upwards or downwards. CTRANS' results for the frac- 
tion of incident flux emerging from each face upwards or 
downwards are compared with those of McKee and Cox in Table 
3.2.1 for each face for upward-emerging and downward-emerging 
photons. Results for all three solar zenith angles are com- 
puted simultaneously so that there will be some correlation 
between results for different suns; suns at 30° and 60° 
correlate most strongly. 

From the results of Table 3.2.1, the conservation of 
flux can be checked: for solar zenith angles, 0°, 30°, and 

60°, the total emerging flux is .946 +_ .12, .98 + .14, and 
1.29 +_ .25. In each case, the difference from 1.0 is not 
statistically significant. 

Overall, the results for individual faces tend not to 
be as accurate as might be desired. For a number of faces, 
the difference between CTRANS' results and those of McKee and 
Cox is large compared with the computed standard error and 
in most of these instances, CTRANS' results are low. This 
may be caused indirectly by photons escaping from the cloud 
before they have added their full proper contributions to 
the flux. Two possibilities exist for correcting this: 
keep the photons in the cloud and accumulate the probability 
that they would have escaped; or run many more photons. The 
first method precludes interaction with the atmosphere sur- 
rounding the cloud and with other clouds and is therefore 
generally undesirable for CTRANS. The second method will 
only be feasible for optically thin clouds. One other possi- 
bility exists which has not yet been fully explored: distri- 

bute the photons non-uniformly over the cloud face initially. 
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(COX-McKEE) 







This is tricky, however, as great care must be employed to 
avoid introducing bias in the final result. Further, a 
suitable non-uniform distribution would be expected to 
depend upon the solar zenith angle, thus limiting runs to 
a single sun. The logical extension of this concept would 
be to attempt to distribute photons both in position and 
direction according to their importance in the final results. 
Initial effort along this line were not notably successful 
and have not been pursued. 


3.3 FLUX COMPUTATIONS FOR A FINITE CLOUD COMPUTED BY A 

FORWARD TRACKING PROGRAM 

Our experience in applying the backward tracking Monte 
Carlo program to the evaluation of fluxes emanating through 
the faces of finite clouds demonstrated that there is con- 
siderable difficulty in establishing relevant importance 
functions which provide unbiased solutions with low variance. 
A more natural procedure for evaluating the fluxes through 
the faces of a single cloud is to use the forward tracking 
procedure of Cox and McKee. 

In order to provide a means of testing flux results 
generated by the backwards code we modified an existing for- 
ward tracking code for a plane parallel atmosphere using sub- 
routines developed for the backtracking code. The code was 
tested by regenerating some of the results of Cox and McKee. 

The flux evaluation code, FLX, provides a capability 
in simulation capacity that exceeds that of the Cox and McKee 
code. The salient features of the simulation are as follows: 
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1. Only one box-shaped cloud can be treated. Its 
interactions with the surrounding atmosphere, 
other clouds and the ground are neglected. The 
cloud must have its upper and lower surfaces 
parallel to the ground. Its shape need not be 
cubic . 

2. The cloud can have arbitrary azimuthal orienta- 
tion, i.e., it is not necessary for there to be 
a pair of faces parallel to the solar zenith 
plane . 

3. Vertical inhomogeneities within the cloud are 
allowed. Up to 100 vertical layers may be 
defined with up to 5 scattering species. 

4. Polarization is treated using the Mueller 
scattering algebra. 

In brief, the photon tracking and sampling algorithm 
is as follows: 

A photon enters a face of the cloud in a random 
location. The face of entry is selected on the basis of 
the area of the face projected perpendicular to the solar 
direction. At each scattering we compute the face towards 
which the photon is travelling using the subroutine CLOUT 
described elsewhere in this report and evaluate the optical 
thickness x^ to the face. We then sample the quantity 


S • = w. I . e 

x 11 


-T - 
1 
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where w^ = weight of the photon at the ith scattering 

and 1^ = intensity component of Stokes vector at ith 

scattering as the contribution to the 
appropriate flux. 

-t - 

The photon weight is then reduced by a fraction (1 - e 1 ) , 
i . e . , 


w 


i+1 


w A (1 - e Ti ) 


and the photon is forced to make its next scattering prior 
to exit from the cloud. If t. is the actual optical distance 

to the next scattering then t. is distributed as e but with 

the constraint that t. < t.. s 

1 s 1 

The sampled fluxes from the sides of the cloud are 
discriminated in terms of up-welling and down-welling photons. 
Obviously, only upwelling photons contribute to the flux from 
the upper face of the cloud and only downwelling photons con- 
tribute to the flux from the lower cloud face. By treating 
the point of entry to the cloud as a scattering point we are 
able to properly evaluate the flux contribution of attenuated 
but unscattered solar radiation contributing to the flux. 

Each photon is tracked until its weight is reduced to 
an arbitrary threshold. Typically we have set the threshold 
at 10 ^ . A typical simulation would involve 1000 to 10000 
incident photons and standard errors in the fluxes range 
from 1 to 5 percent. 
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Table 3.3.1 presents a comparison of the flux fractions 
emanating from the various faces of the cloud. Contributions 
from unscattered light have been added into the diffuse con- 
tribution in computing the total downwelling flux. TINF and 
BINF refer to corresponding results for a cloud layer of 
infinite horizontal extent. Apart from these differences. 
Table 3.3.1 is similar in format to Table 3.2.1 and permits 
comparisons among the results from CTRANS , FLX, and McKee 
and Cox. 


3.4 GROUND PLANE REFLECTIONS 

The basic mechanism for handling Lambertian reflection 
has been tested implicitly in the Rayleigh atmosphere tests 
reported in Section 3.1. The capacity of CTRANS to handle 
inhomogeneous Lambertian ground patterns may be illustrated 
by a test example but specifically comparable independent 
calculations do not appear to be available. 

As a suitable test example, the basic ground patterns 
were chosen as half planes of differing albedo with boundaries 
parallel to the y-axis and x-axis intersections at a number 
of points along the axis. The receiver was above the top of 
the atmosphere (704 km) centered at (xj^, Y^) = (0,0) and 
looking towards the nadir. The receiver field of view was 
85.2 yrad yielding a 60 m diameter footprint on the ground. 

The atmosphere was composed of Rayleigh and haze C aerosol 
scatterers with total optical thicknesses of .098 and .19643, 
respectively. The wavelength employed was X = 0.55 ym . The 
albedo pairs employed were (.6, 0), (.5, .1), (.4, .2), 

(.3, .3) with boundaries intersecting the x-axis at 
X = ^2.0, _+l . 0 , ^0.5, +0.25, +0.05, 0.0 km. In addition to 
the half-planes, results for uniform Lambertian ground planes 
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with albedos .6, .5, .4, .3, .2, .1, and 0.0 were obtained 
All results were, of course, obtained simultaneously. These 
results are shown in Figure 3.4.1. Since all computed values 
depend upon the same set of photon tracks, relative values 
are computed more accurately than absolute values. This 
feature is of great importance when CTRANS is used, for 
example, to determine the distinguishability of surface 
features as seen from a satellite. 

The ground plane reflection characteristics may be 
modified in yet another way: the reflectance type may be 

changed from Lambertian to rough Fresnel, useful for sea 
surface modeling. Several independent calculations exist 
against which CTRANS' results may be validated. We chose 
the following parameters to specify a suitable test case: 
solar zenith angle, 0^=57°; wind speed=5 m/sec (yields a 
wave slope variance a 2 =0 . 0286); wavelength , A=0.70 pm; 
Rayleigh optical thickness x^=.037; aerosol optical thick- 
ness t m = 0.214; the field of view of the receiver was taken 
to be zero. These parameter values were chosen so as to 
be able to compare our results with those computed inde- 
pendently by RADTRAN. These model parameters are also 
similar to (but not identical with) those used by Kattawar 
and Plass . ^ Figure 3.4.2 presents the results of 

CTRANS as well as those of RADTRAN and Kattawar and Plass. 

As may be seen, all three calculations are in agreement 
to within 1.4 times the computed standard error. 
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FIGURE 3.4.2. 



SECTION 4.0 

SUBROUTINES: DESCRIPTIVE OUTLINES 


In its present version, CTRANS is composed of a MAIN 
control routine and 34 supporting subroutines. In the follow- 
ing, we outline the main functions of each of these. 

MAIN 


MAJOR FUNCTION: Main control routine 


INPUTS: 

1 . 

Initial random number 


2. 

IND, 1ST, I SKY, INFLUX, WAV 

OUTPUTS: 

1 . 

Contribution maps 


2. 

Timing information, number of photons tracked 

ACTIONS: 

1 . 

Set up constants, initiate random number 
sequence, clear storage. 


2. 

Read major control card: IND, 1ST,... 

(STOP if IND=0) 


3. 

Call INPUT - Input model specification, 
write out descriptive text 


4. 

Loop over photon tracks 

4.1 Initialize photon at receiver 


4.2 Sample direct sun contribution (use 
SEESOL) 

4.3 Loop over allowed number of scatterings 

4.3.1 Determine distance to scattering, 
propagate photon along path (TRACK) 

4.3.2 Determine scattering species 

4.3.3 Sample hypothetical photon from 
each sun (SAMPLE) 
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4.3.4 


Determine scattering angle and 
direction (ANGLEM and JIM; or 
GNDREF) 

4.3.5 Update cumulative scattering 
matrix following a real scat- 
tering (RENEW) 

4.4 Clump, collect statistical measures 

4.5 Check time remaining (terminate if 
insufficient) 

5. Write contribution maps (if finite receiver 
field of view) 

6. Normalize samples 

7. Output results (OUTPUT) 

8. Output statistics (TAB) 

9. Return to Step 2. 
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INPUT 


MAJOR FUNCTION: Read in model specification parameters, 

compute and fill reference tables, write 
text describing model and input quantities. 


INPUTS: 1. 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 
9. 


Atmosphere description cards 

Sun specifications 

Size distribution parameters 

Phase matrix elements (cards or tape) 

Atmospheric layer heights, temperatures 

Atmospheric layer heights, component densities 

Cloud input specifications (CLINP) 

Receiver specifications 

Ground reflectance type map specifications 


OUTPUTS: 1. 

2 . 

3. 

4. 

5. 


6 . 


Main title 

Simulation parameters 
Receiver description 
Solar array description 
Ground description 

5.1 Reflectance type map 

5.2 Lambertian albedo maps 
Atmosphere description 

6.1 Characteristics of each type of scat- 
terer in entire atmosphere 

6.2 Phase functions 

6.3 Ambient atmospheric densities, cross 
section 

6.4 Average scattering angle 

6.5 Cloud specifications 
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SAMPLE 


MAJOR FUNCTION: Accumulates a sample for hypothetical 

photons arriving from each Sun. Reads 
no cards. Writes no text. 

ACTIONS: 1. Set up photon co-ordinates (if in a cloud) 

(TRANS 1) 

2. Loop over Suns 

2.1 Compute scattering angle and/or direc- 
tion for scattering into Sun 

2.2 Determine scattering plane, normal to 
scattering plane, rotation angle (FIXALL) 

2.3 Multiply rotation into cumulative 
scattering matrix (ROTATE) 

2.4 Obtain matrix elements for scattering 
other than Rayleigh (ELEMTS or GENGND) 

2.5 Multiply scattering matrix into cumula- 
tive scattering matrix (SCATMI , SCATRA 
or GENSCA) 

2.6 Compute optical path length to Sun 
(SEESOL) and, thence, the attenuation 
factor 

2.7 Fold in other applicable weights 

2.8 Accumulate sample into sampling matrix, 

VS AMP 
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SEESOL 


MAJOR FUNCTION: Determines optical path length to the Sun 

Assumes that CLCORD and COORD contain 
position and path direction cosines for 
photons in a cloud or not in a cloud, 
respectively. 


ACTIONS: 


1. Determine distance out of cloud (if in a 
cloud) (CLOUT) 

2. Compute contribution of overlying ambient 
atmospheric layers 

3. Compute contribution of present ambient 
atmospheric layer 

4. Compute contribution of other clouds entered 
along path 

4.1 Check cloud spheres entered along path 
(CLENTR) 

4.2 For each sphere entered, check entries 
in detail; get distance to entry and 
distance to exit for each cloud entered 
(SDIST, CLOUT) 

5. Compute total optical distance to Sun (along 
the direction indicated by COORD and/or 
CLCORD) 
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TRACK 


MAJOR FUNCTIONS: 1. Determines optical path length to 

be traversed 

2. Propagate photon along its current 
direction until chosen optical path 
length is traversed 

3. Update co-ordinate systems at boundary 
(cloud) crossings 

4. Set IOUT=l for photons exiting top of 
atmosphere 

5. Set ISCAT=1 for photons hitting earth, 
compute coordinates of impact 

ACTIONS: 1. Choose optical path length to be traversed 

x = -5,n(P) (P=uniformly distributed random 
number) 

2. Set up boundaries 

2.1 Boundary of present cloud (if in a cloud) 

2.2 If not in a cloud, boundary of any cloud 
along extended path - use CLENTR, TDIST- 
or, atmospheric boundaries if no clouds 
can be entered 

2.3 Determine type of boundary 

3. Successively accumulate increments in distance 
traversed and traversed optical path length 
until either the chosen optical path is ex- 
hausted or until a boundary is encountered. 

4. Boundary encounter 

4.1 If boundary is top of atmosphere, set 
I0UT=1, exit 

4.2 If boundary is ground, set ISCAT=1, com- 
pute coordinates of impact, CG. 

4.3 If boundary is a cloud boundary, update 
coordinates and resume tracking 

5. Update photon position 
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TRANS 

TRANS1 

TRANS2 


FUNCTIONS: These closely related and similar routines 

transform coordinate positions and/or photon 
flight direction cosines into and out of the 
cloud coordinate systems. 

If XC^ = photon position coordinate in the 
cloud i=l,3 

X^ = photon position in fixed coordinate 
system i=l,3 

XC^ = photon flight direction cosines in 
the cloud i=4,6 

X^ = photon flight direction cosines in 
the fixed coordinate system i=4,6 

and XCL^ = C.M coordinates of the cloud in the 



the fixed system 

i = l , 3 

0 

= orientation 

angle 

of the cloud 

xc 

= R(X-XCL) 

i = l , 3 


xc 

= R(X) 

i=4 , 6 


X 

= XCL + R XC 

i = l , 3 


X 

= R XC 

* 

II 

•H 



/ cos0 sin0 

°\ 


R 

= |-sin0 cos0 

0 



\ 0 0 

if 


R 

= transpose of 

R 
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CLOUT ( I , NOUT , CD I ST , EXT) 


FUNCTION: 


PARAMETERS: 


METHOD: 


Determines the point of exit from a cloud and 
the distance along the photon path to this 
point 

I = Cloud index 

NOUT = 1 if exit point is found 

CDIST(I)= distance to exit point 

EXT(I , J) , J=1 , 2 , 3 = exit point coordinates for 

cloud I 


Let C. , C . ,C, 
l * 3 * k 


WA 


v i’V v k 


W A k 


exit coordinates in cloud 
system 

present photon coordinates in 
cloud system 

direction cosines of photon 
flight path in cloud system 
cloud half edge lengths 
all where i,J,K may be a per- 
mutation of X,Y,Z. 


= V. *A. 

l 1 l 


C. = X.+(C.-X.)V./V. 
3 3 113 1 

X, +(C.-X. )V v /V. 
k v 1 1 ' k' 1 


C, = 


if C- and C r - are within the cloud 
3 k 

boundaries then the exit point is 

(C.,C.,C-,_) i if not, cycle through 
1 3 K 

values of (i , j , k) 

Compute distance to exit point 
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CLINP 


FUNCTIONS: 


ACTIONS: 


Reads in cloud input parameters, transforms 
these and computes auxilliary descriptive 
parameters . 


1 . 

2 . 

3. 


4. 


Read number of clouds NCLOUD 
Return if NCLOUD=0 
For each cloud read 


HC , AC , BL , CC , XCL , YCL , THETAC , INC , CLDEN 


where 

HC=Altitude of base of Cloud (in km.) 

AC , BC , CC=Linear dimension of cloud (km.) 
XCL, YCL =X,Y coordinates of center of 

mass of the cloud in the fixed 
coordinate system 

THETAC =Orientation angle of the cloud 
measured clockwise looking up 
(in degrees) 

INC =Indicator of scatterer type 
composing cloud 

CLDEN =Density of scatterers of type 
INC (cm -3 ). 

Computes sines and cosines of THETAC, 
cloud edge half-lengths, radii of cloud 
circumscribing spheres and puts all 
information into the descriptive vector 


CLOUD (10,10) : 

CLOUD(I,J) : I=Cloud pointer 
J=l,2,3 : C.M. coordinates of the cloud 

J=4,5,6 : Half dimensions of each side 

of the cloud 


J= 7 , 8 
J=9 


J=10 


: cos (THETAC) , sin (THETAC) 

: Radius, R, of the circumscribing 
sphere 
: R 2 
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CLENTR 


FUNCTIONS: 


ACTIONS: 


Checks whether the photon can penetrate any 
of the cloud circumscribing spheres and orders 
the sequence of entry. 

1. Computes square of distance to a cloud, 
length between point of closest approach 
and present photon position, and square 
of distance of closest approach along the 
extended path. 

2. Orders entries by entry sequence 
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SDIST 


FUNCTIONS: 

ACTIONS: 


TDIST 

FUNCTION: 

ACTION: 


Tests which plane of the cloud photon entered, 
calculates distance to intersection, and orders 
them by entry sequence 

1. For each uneliminated cloud computes point 
of entry and distance to that point (PLANE) 

2. Orders by entry sequence (length from 
present position to entry point. 


Calculates distance to the cloud photon can 
enter in TRACK 

Uses PLANE to compute distance to the first 
cloud entry 
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PLANE ( I , INDC , ENT , NPLANE , DIST) 


FUNCTION: Calculates the coordinates of the point of 

entry of a photon into a cloud and determines 
the distance to this point from the photon’s 
present position 

PARAMETERS: I = Cloud index 

INDC = Parameter controlling the transforma- 

tion into/out of cloud coordinate 
system 

ENT(I,J) J=l,2,3 = Coordinates of point of 

entry into cloud I 
NPLANE = 1 if entry can occur 
= 0 otherwise 

DIST(I) = Distance to point of entry into 
cloud I 

METHOD: Let C^,Cj,C^ = Entry coordinates in cloud system 

X-,X-,X V = Present photon position coordi- 
1 1 K 

nates in cloud system 

V-,V-,V v = Direction cosines of photon flight 
1 j K 

path in cloud system 
A.,A.,A V = Cloud half-edge lengths 

all where i,j,k are a permutation 
of X , Y , Z . 


1. 

If 

X. 1 <A - 
1 r— i 

try another permutation 

2. 

If 

X - <A. 

C 

. =-A. 



l l 


l l 



X. >A. 

C 

. =A. 



l l 


l l 

3. 

If 

(X. -C. )A. 

v l i' l 

>0 1 



or 

v.=o. 


| Try another permutation 



1 



4. 

c. 

3 

= x. + 

3 

tv s 

/V i ) (C i -x i ) 


c k 

■ A * 

(A 

/v i ) (C i -x i ) 
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FIXALL (UX, UY, UZ , CX , CY, CZ ,UN ,VN ,WN, ANGLE , A , B, C , ) 


FUNCTIONS: 


PARAMETERS: 


ACTIONS: 


Determines normal to the scattering plane 
and the angle through which the Stokes 
vector must be rotated 


U = UX , UY , UZ 

C = CX , CY , CZ 

N = UN , VN , 1VN 

ANGLE 
M = A, B , C 


Direction cosines of post- 
scattering photon 
Direction cosines of pre- 
scattering photon 
Direction cosines of normal 
to previous scattering plane 
Necessary angle of rotation 
Direction cosines of normal 
to new scattering plane 


1. Checks (UX,UY,UZ) and CX,CY,CZ) for 
parallelism if parallel, return with 
(A , B , C) = (UN , VN ,WN) , ANGLE=0 . 

2. Computes direction cosines of normal to 
scattering plane, (A,B,C) via NORMAL 

3. Computes dot product N*M and thence 
the rotation angle magnitude 

4. Sets ANGLE=- ANGLE IF C-(MxN)<0. 
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NORMAL (CXI, CY1, CZ1.CX2, CY2.CZ 2, CX,CY,CZ) 


FUNCTION: 


PARAMETERS: 


ACTION: 


Computes direction cosines of unit vector 
normal to the plane defined by two vectors 


VI = (CXI , CY1 , CZ1) Direction cosines of one 

vector defining plane 

V2 = (CX2 , CY2 , CZ2) Direction cosines of other 

vector defining planes 

N = (CX,CY,CZ) Direction cosines of NORMAL 

1. Computes N = VI x V2 

2. Normalizes N to unit magnitude 

VlxV2 

i.e. N = 

VlxV2 
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JIM (COUT ,CIN,CP,SP,CB,SB) 


FUNCTION: Computes photon direction cosines after 

scattering given the direction cosines 
prior to scattering and the polar and 
azimuthal scattering angles 


PARAMETERS: 

COUT ( J) 

J=l,3 
J=4 , 6 

Photon position coordinates 
Post-scattering direction 
cosines 


CIN(J) 

J=1 , 3 
J=4 , 6 

Photon position coordinates 
Pre- scattering direction 
cosines 


cp,sp 

Cosine 

tering 

and sine of azimuthal scat- 
angle 


cb,sb 

Cosine 

angle 

and sin of polar scattering 

ACTION: 

Let 0 O 
then 

= polar 

coordinates of CIN(J), J=4,6 


COUT (4) = CB*CIN(4)-SB*SP*SIN(4>) + SB*CP*CIN(6)* 
COSO) 

COUT (5) = CB*CIN(5) + SB*SP*COSO)+SB*CP*CIN(6)* 
SIN O) 

COUT (6) = CB*CIN(6)-SB*CP*SIN(0) 
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XSECTM 

XS2 

FUNCTION: 


METHOD: 


These essentially similar routines perform 
table look-up and interpolation of the phase 
function given a value for the scattering 
angle, B. XSECTM operates only on S^ and 
XS2 operates on general phase matrix 
elements . 

Three point Newton's method - see for example 
R.W. Hamming, Numerical Methods for Scientists 
and Engineers , McGraw-Hill, New York, 1962; 
p 99 . 
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ANGLEM 


FUNCTION: 


PARAMETERS : 


METHOD: 


Given a rectangularly distributed random 
number, and the scattering species index, 
chooses a scattering angle by table look-up 
Uses the scattering angle distribution func 
tion formed from for the indicated 

species . 

Where JK points to the scattering species 

X1(JK,100) = Table of scattering angles for 
species JK 

Y1(JK,100) = Probability distribution func- 
tion for the angles listed by 
XI computed from S-^. 

2 or 3-point Newton's method 
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ELEMTS (BETA, S, K) 


FUNCTION: 


PARAMETERS : 


ACTIONS: 


Computes all the scattering matrix elements 
(for an atmospheric scattering) for a 
specified scattering angle BETA. 

BETA = Scattering angle 

S(I) = Result value of I 11 * 1 matrix element 

S(l) = s n 
S (2) = S 12 
s (3) . s 33 
S (4) = S 34 
S(5) = S 22 
S(6) = S 44 

Use XSECTM to obtain S(l) for hypothetical 
scattering (S(l)=l for real scattering). 

Use XS2 for remaining matrix elements. 
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RENEW (COUT , BETA , CB) 


FUNCTION: 


PARAMETERS: 


ACTIONS: 


Updates the cumulative scattering matrix 

following a real scattering 

COUT = Six-vector describing post- scat tering 
photon 

BETA = Scattering angle 

CB = cos (BETA) 

1. Determine normal to new scattering plane 
and the angle between it and the normal 
to the previous scattering plane (or 
reference direction) using FIXALL 

1.1 Update scattering plane normal 

2. Post- multiply a rotation into the 
cumulative scattering matrix 

3. Obtain scattering matrix elements and 
accumulate a scattering by post-multipli- 
cation 

3.1 For ground scattering obtain 
matrix elements from LAMBRT (Lam- 
bertian scattering) or from GENGND 
(rough Fresnel scattering) and use 
GENSCA to accumulate the scattering. 

3.2 For atmospheric scattering, obtain 
Mie or generalized scattering matrix 
elements from ELEMTS. Accumulate 
the scattering using SCATRA (Rayleigh 
scattering-matrix elements computed) , 
GENSCA (general scattering) , or SCATMI 
(Mie scattering) . 

4. Update CIN to the values contained in COUT. 
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FUNCT (A1,B1,C1,D1,IDIST,D) 


FUNCTION: Computes the value of the particle diameter 

probability density function for a given 
particle diameter. 


PARAMETERS: A1,B1,C1 

= distribution function 

parameters 

D1 

= Normalization factor 


IDIST 

Specifies type of distribution 


= 1 

Junge 



= 2 

Deirmendj ian 



= 3 

Log-normal 



= 4 

Haze-C (Junge with 

pedestal) 

D 

= Particle diameter 



USE: FUNCT is called as a function 


ACTIONS: 

IDIST=1 

FUNCT = D1*D‘ B1 

D>A1 



= 0. D<A1 



IDIST= 2 

FUNCT = D1*D A1 

exp (-B1*D^^) 


IDIST= 3 

FUNCT = 0. 

D<C1 


D- Cl eXp ^ ^ A1K2,n( -Bl-Cl^ * 

IDIST=4 FUNCT = 0 D<C1 

= D1 C1<D<A1 

at BI 

= D1 (g— ) D>A1 
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GENSCA (SSIN , SOUT , SK , K) 
SC ATM I (SSIN, SOUT, S,K) 
SCATRA (SSIN, SOUT, CB,K) 


FUNCTIONS: 


PARAMETERS : 


ACTIONS: 


These .similar routines accumulate a new 
scattering by post-multiplying the scatter- 
ing phase matrix into the cumulative scat- 
tering matrix for real or hypothetical 
(sampling) scattering. 

SSIN = 16-element vector containing the 
elements of the 4x4 cumulative 
scattering matrix arrayed column- 
wise . 

SOUT = 16-element vector of elements of 

cumulative scattering matrix after 
scattering. 

S = 6-element vector of matrix element 
values 

K =1 Compute only first column of SOUT 
(i.e., first 4 elements) for 
hypothetical (sampling) scattering. 
All elements are not needed for 
unpolarized incident light. 

K =2 Compute all elements of SOUT. 

CB = Cosine of scattering angle used by 

SCATRA to compute scattering matrix 
elements . 

If S , is the output matrix, Sj is the 
input matrix 
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SPECIALIZATION: SCATMI Assumes that S^=S^ and S-^=S 22 

Assumes that S^=l if K=2 

SCATRA Computes S n =| (1+ (CB) 2 )=S 22 

S 12 =f ((CB) 2 -!) 


S 33 = ! CB=S 44 

s 34=°. 

Takes S^ = l anc ^ normalizes other 

matrix elements by 
dividing by ^-(1+(CB) 2 ) 


if K= 2 . 
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ROTATE (SSIN, SOUT,PHI) 


FUNCTION: 


PARAMETERS: 


ACTIONS: 


Folds a rotation into the cumulative scatter- 
ing matrix by post-multiplication. 


SSIN = 16-element vector containing column- 
wise ordered elements of the cumulative 
scattering matrix 

SOUT = 16-element vector containing elements 
of the cumulative scattering matrix 
after rotation. 

PHI = Rotation angle 


S q is the cumulative scattering matrix 
= after rotation 

Sj is the cumulative scattering matrix 

= before rotation 

4> = rotation angle 


S 


o 


r - 



0 

cos 2cf) 
-sin2<j) 
0 


0 

sin2cj> 
cos 2 4> 
0 
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REFLOC (CG, IRCOD) 


FUNCTION: 


PARAMETERS: 


ACTIONS: 


Returns a value of IRCOD indicating the 
reflectance type at a point indicated by 
CG corresponding to the map choice fixed 
by INPUT. 

CG(1),CG(2) = X-Y location of point of impact 
of photon on the ground 
IRCOD = Reflectance type code 

= 0 Lambertian 
= 1 Rough Fresnel 

1. If ground is uniform plane (i.e., 

IRPAT=0) IRCOD=IRBGD 

2. If map is half-planes |IRPAT|=1 
then 

for IRPAT=+1 and X>XBND IRCOD=IRSPOT 
for IRPAT=- 1 and X<XBND IRC0D=IRSP0T 
with 

IRCOD=IRBGD otherwise 

3. If map consists of rectangles differing 
from background (IRPAT=2) 

then 

IRCOD= IRSPOT if | X- CNTRX | <EDGEX/ 2 
and | Y-CNTRY| <EDGEY/2 
IRCOD=IRBGD otherwise 
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LAMBRT (COUT , S , K) 


FUNCTION: 


PARAMETERS: 


ACTIONS: 


Constructs matrix elements appropriate to 
Lambertian scattering 

COUT = Descriptor vector for scattered photon 
C0UT(6)= zenith cosine 
S = 6-vector of matrix elements 

K =1 Indicates hypothetical scattering 

K =2 Indicates real scattering 

S(l) = 2*COUT(6) K= 2 

SCI) = 4* COUT (6) K= 1 

S(J) = 0 J= 2 , 6 
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GNDREF (COUT , CIN , XTH , XPH , IRCOD) 


FUNCTION: 


PARAMETERS: 


ACTIONS: 


Generates a reflection direction from the 
ground 

COUT = Post-scattering photon descriptor 
CIN = Pre- scattering photon descriptor 
XTH, XPH - Uniformly distributed random 
numbers 

IRCOD = Reflectance type indicator 
= 0 Lambertian 
= 1 Rough Fresnel 

1. For Lambertian reflection 
<p = 2tt*XPH 

Cos (0 ) = XTH 

COUT (4) = SIN0 COScj) 

COUT (5) = SIN0 SINc p 

COUT (6) = COS0 

COUT ( J) = CIN (J) , J=1 , 3 

2. For rough Fresnel reflection 

2.1 Find slope X,Y components from 
Gaussian distribution through GAUSS 
(standard deviation=WS IG, mean=0.) 

2.2 Compute normal to slope (direction 
cosines labeled WNG(I), 1=1,3) 

2.3 Check to see if slope is visible 

to incident photon (visible if 
3 

G = £ CIN (J) *WNG ( J) <0 . ) 

J=1 

2.4 COUT (J+3) = - 2G*WNG (J) + CIN(J+3) forms 

the direction cosine of 
reflection direction 
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GENGND (COUT , CIN 
FUNCTION: 

PARAMETERS: 


ACTIONS: 


, IRCOD , S', K) 

Computes Fresnel matrix elements multiplied 
by statistical factors 

COUT, CIN = Scattered and unscattered photon 
descriptors, respectively 
IRCOD = Not used 

S = 6-vector of matrix elements 

K =1 Hypothetical (sampling) scat- 

tering 

= 2 Real scattering 


Let V q and be the (unit) photon propagation 

vectors following and prior to reflection, 
respectively and n be the normal to the 
surface facet. 

1. Compute cosine of angle of incidence: 


cos X 




2. Compute factor to account for non-vertical 
incidence of photon (modification to 
probability of encounter for a chosen 
facet slope) 

FACTR = cos X/[(n*Z) (-Z*Vj)] 

3. Compute cosine of angle of refraction 


cos 9^ = Jl - (1-cos^X) f =Index of 

L f refraction 


4. Compute Fresnel factors 
f cosX - COS0^ 


R 


11 


f cosX + cos6. 


R i ■ 


cosX - f COS0 


cosX + f COS0 
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Compute matrix elements 

CO 

n 

( R b 

+ Rp/2 = S 2 

S 1 2 

( R n 

- Rj)/2 

S 33 

S 44 " 

R 1 R 11 

S 34 = 

0 



6. If K=l, compute probability of encountering 

proper slope ; FEXP 

7. Fold all auxilliary factors into matrix 
elements . 
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TAB (VSAMP2 , VS AMP 3 , XMC , XMX , IND , IFL) 


FUNCTION: 


PARAMETERS : 


ACTIONS: 


Compute and display mean, variance, standard 
deviation, and standard error for each sun, 
each Lambertian albedo pattern, and each 
Stokes’ parameter (I,Q,U,V) for single 
scattering, multiple scattering and total. 
Compute correlation coefficients amongst 
suns . 


VS AMP 2 


VSAMP3 

XMC 

XMX 


IND 

IFL 


Sample matrix VSAMP2 (IQUV, ISUN , IS , 1 , IALB) 
where 

IQUV Points to Stokes' parameters 
ISUN Designates the sun 

IS Designates direct sun contribution, 
single scattering, and higher order 
scattering . 

Sample matrix of squares of clumps 
Number of clumps 

Factor used to normalize results to 
the number of photons actually tracked 
(when this varies from NMAX) . 

5 or 6 display results only for zero 
albedo 

Indicator designating whether answers 
are for flux or intensity 


1. Compute normalization factors 
f=(solar const. )/4tt) 

2. Compute and write means and standard 
deviations 

std. dev. = [ VS AMP3*f 2 *XMC- (MEAN) 2 ] 

2 2 

3. Compute variance = (std. dev.) = a 
and standard error = Vc^/XMC 

4. Compute correlation coefficient matrix 
for results for each solar direction. 
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OUTPUT (VSAMP , IND , FACT) 


FUNCTION: 


PARAMETERS : 


ACTIONS: 


Normalize results, write out intensity or 
flux results labeled by Lambertian albedo 
pattern. If appropriate, compute effective 
albedo evaluated from regression on uniform 
albedo patterns and modulation transfer 
functions . 

VSAMP = Sample matrix 
IND = Not used 

FACT = Normalization adjustment factor 
(normally FACT=.25) 

1. Compute normalization factor 
f=(solar const .) *FACT/tt . If receiver 
is the face of a cloud, normalize 

to unit flux incident on the cloud 
(Cox-McKee normalization). 

2. Write out results (normalized) for 
intensity or flux. For intensity, 
results are direct sun contribution, 
single scattering, higher order scat- 
tering, and total intensity. For flux, 
results are direct solar contribution, 
upward flux, downward flux, total flux. 

3. If the Lambertian albedo patterns con- 

tain uniform planes with A=0. and at 
least two others (and if the reflectance 
type map is uniformly Lambertian) perform 
a regression analysis: find AL,BL in the 

expression 

[ I (A) - I (0) ] /A = AL + BL* [I (A) - I (0) ] 

where A refers to the albedo and 1(A) 
is the result for albedo=A. Use AREG. 
Also computed are the effective albedos 
for other albedo patterns. 
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4. If appropriate (IMTF=1, set by ALBEDO) 
compute percent contrast from pairs of 
results both for intensity and effective 
albedo . 
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ALBEDO ( IALB , ALBPT) 


FUNCTION: 


ACTIONS: 


Write titles or compute the albedo at the 
point specified by CG. 

. IALB = -J Write out description (50 char- 
acter maximum width) of albedo 
pattern J. 

IALB = +J Compute albedo at the point 

X=CG(1), Y=CG(2) and return this 
value as ALBPT. 

. If this pattern is for computing MTF's, set 
IMTF=1 . 
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AREG (NALB , AL , BL) 


FUNCTION: 


PARAMETERS : 


ACTIONS: 


Perform albedo regression analysis assuming 
FTOT (NALB) is the value corresponding to a 
uniform zero albedo plane. 


NALB = Total number of albedo patterns 
NUNIF = Number of uniform plane albedo 
patterns 

AL,BL = Regression coefficients 


1. Use LINFIT to compute albedo regression 
coefficients AL,BL 

Fo?m : S 5 i0n [ I (A) - I (A=0 . ) ] /A=AL+BL* [ I (A) - I (A=0 . ) ] 

2. Compute effective albedo (AEFF (I ) , 1=1 ,NALB- 1) : 


AEFF(I) 


I (J) -I (NALB) 

AL + BL* (I (J) - I (NALB) ) 


AEFF (NALB) = 0. 
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LINFIT (N,Y,X,A,B) 


FUNCTION: 


METHOD : 


Performs lines 
parameters A,E 

Y=A+B*X where 
length N. 


n E 
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SECTION 5.0 
INPUT. PARAMETERS 


The input parameters are listed below, together with a 
concise description of their use, meaning, and default values 
(if applicable) . The order of appearance is that of a proper 
data deck. 

Initial random number (may be blank, default value = 65549) 
(17) 

IND , I ST , I SKY , INFLUX , WAV 
(411, F8. 5) 

IND = Calculation extent indicator 
= 0 Stop calculation 

= 2 Input new clouds and receivers only 
= 4 Normal start 

= 5 Omit polarization, set scattering matrix 
S i j=0 for i , j 1 

= 6 Calculate only for completely absorbing ground 
= 7 Completely absorbing ground + no polarization 
1ST = New atmosphere indicator 
= 1 Read in new atmosphere 
= 0 Use one already read in 
ISKY = Sun indicator 

= 0 User specified suns (maximum of 5) 

= 1 Suns at 9 = 0? ,30? ,60? ;<p = 0 

= 2 Suns at 9 = 0? ,30? ,60? ;43 = tt/4 

= 3 Suns at 9 = 0? , 30? , 60? ; 4> = tt/2 

INFLUX = Indicator specifying whether flux or intensity is 
to be calculated. 

= 0 Calculate flux at receiver 
= 1 Calculate intensity at receiver 
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WAV = Wavelength of radiation (in micrometers) . 

INDIC (5) , IDIST (5) , NLAYER , NSOL , NMAX ,NSCA , NALB ,IUNIT (5) ,NT0T(5) 
(1011,213,16,712,513) 

INDIC(k) = Indicator for type of k*"* 1 species. 

> 5 If species is not present 

= 0 For Rayleigh (matrix elements computed) 

= 1 For monodisperse Mie (matrix elements 

read from binary tape) 

= 2 For polydisperse Mie (matrix elements 

computed while reading from binary tape) 

= 3 For experimental data (matrix elements 

read from cards) 

= 4 For ozone (absorption only) 

IDIST(k) = Particle diameter distribution indicator 

"th 

for k species. 

= 1 For Junge distribution 
= 2 For Deirmendj ian' s Haze model 
= 3 For log-normal distribution 
= 4 For Haze C (Junge with pedestal) 

NLAYER = Number of plane parallel layers in ambient 
atmosphere . 

NSOL = Number of solar zenith angles to be con- 

sidered . 

NMAX = Number of photons to be tracked 

NSCA = Maximum number of scatterings allowed per 

photon track 

NALB = Number of Lambertian albedo maps 

IUNIT(k) = Indicator of (tape) unit from which matrix 
elements of k L species is to be read. 

NTOT (k) = Number of diameter steps to be read for 

■f-L 

the k (polydisperse) scattering species. 
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INDICA (5) 

(511) 

INDICA(k) = Indicator similar to INDIC(k) but applicable 
to ambient atmosphere (not clouds) . 

Note - Any species appearing in INDICA must 
also appear in INDIC. 

SCALE (5) 

(5E11.6) 

SCALE^k) = 0. Use cross section as otherwise read in 

and/or computed for species k. 

>0. Scale cross section of ambient atmosphere 
component k so that the total vertical 
optical thickness of component k=SCALE(k). 

AR,BR 

(2F8.6) AR and BR specify the Rayleigh cross section 

a R ; 

_ o o _ T)n 

o r =ARx10 xX , where X is the wavelength. 

ADEF , EXPFAC 
(2F8. 5) 

ADEF = Ozone cross section 

a n =ADEF/2.687xl0 19 
3 

EXPFAC = Scale height adjustment parameter for 

molecular species. If H is the effective 
scale height computed from densities as 
read from cards, the adjusted scale height 
H'=H/ (1-H*EXPFAC) . 
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INDOZ , INDCOZ 

( 211 ) 

INDOZ = Ozone input units specification 

= 0 Ozone densities cm ^ are read in 
f 0 Ozone partial pressures are read in 
= 1 Ozone densities computed from partial 
pressures are temperature corrected. 
INDCOZ Inactive parameter, no function ascribed 
for it in this version. 

SOLC 

(F10.6) 

SOLC = Solar constant assumed for normalization 

= 0. (or blank) on input, a value of 1.0 
is adopted. 

f 0. Use the specified value of SOLC. 

If ISKY=0, solar directions must be read in as 
follows : 

THETA (1) , PHI (1) 

(F7 . 4 , F7 . 2) 

THETA(NSOL) ,PHI (NSOL) 

THETA(I) = Solar zenith angle for I** 1 sun (in degrees) 

4 . 1 - 

PHI (I) = Solar azimuth from I sun (in degrees) 

NK(1) , STEP (1) ,A1 (1) ,B1(1) ,C1(1) , ITAPE(l) 

(13 , 4F8 . 3 , II) 


89 



[When the K L species' matrix elements are to be read from 
cards, they appear as follows:] 

NK(K) , STEP(K) , . . . 

CSECT(K) , ASECT (K) 

(2E11.4) 

XA,S11(1,K) ,S12(1,K) ,S22(1,K) ,S33(1,K) ,S34(1,K) ,S44(1,K) 
(F6 . 1 , 6E12 . 5) 

XA , $1 1 (NK, K) , S12 (NK, K) , . . . 


NK(KMAX) , STEP (KMAX) ,A1 (KMAX) , . . . (Where KMAX is the maximum 

total number of atmospheric 
species present) 

NK(K) = Number of angles for which the matrix elements 

1 " 

of the K species are to be read from tape 
or cards . 

STEP(K) = Step size in x (x=tt D/A) for polydispersions 
(K tB species) 

= Value of x for monodispersions (K L species) 

A1 (J) , B1 (J) , Cl (J) are the parameters entering into the 

particle size distribution function 
for the scatterer type. 

For IDIST(J) = 1, the diameter distribution is 
F = D1 D' B1 D > A1 

where D1=(B1-1) A1 (B1_1 ^ 

For IDIST (J) =2 

AT H 

F = D1 D A1 exp ( - BB1*D U1 ) 


where D1 = 


Cl B d A1 + 1 ’ /C1 


BB1 = Bl*2 
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For IDIST (J) =3 


F = 


D1 

D- Cl 


exp {- [ A1 log ( 


D- Cl 
Bl-Cl 


2 

)] >, 


D>C1 


= 0. D < Cl 

where D1 = A1//jt 


For IDIST (J) =4 

F = 0. 

F = D1 

Ai B1 

F = D1 (gi) 
where D1 


D<C1 

Cl < D < Al 
Al < D 
Bl-1 

( B 1 - 1 ) ( Al - C 1 ) + A 1 


ITAPE(J) = 0 Use NTOT(J) and STEP(J) as read from cards. 

= 1 Read NTOT(J) and STEP(J) from Mie tape. 


XH (1 ) , TEMP (1 ) 

(F3.0,2X,F3.0) 

XH (NLAYER) , TEMP (NLAYER) 

t ll 

XH(J) = Layer top height of j ambient atmospheric 
layer (in km.) 

TEMP(J) = Temperature of layer (in °K) 

XH (1 ) , DENS (1,1) 

(F3.0,2X,E9.3) 

XH (NLAYER) , DENS (1 , NLAYER) 
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XH (1) , DENS (2,1) 


XH(NLAYER) , DENS (KMAX , NLAYER) 

XH(J) = Layer top height of J ambient atmospheric 

layer (in km.) 

DENS(K,J) = Density of K*'* 1 species in layer (cm 

NCLOUD 

(13) 

NCLOUD = Number of clouds to be considered. 

*** The following are read in if NCLOUD>0 

HC (1) , AC (1) , BC (1) ,CC(1) , XCL (1) ,YCL(1) , THETAC (1) ,INC(1) ,CLDEN(1) 
(F8.5,SF8.3,F8.5,I3,E9.3) 


HC (NCLOUD) , AC (NCLOUD) , . . . 
HC (I) 

AC (I) , BC ( I ) , CC ( I ) = 
THETAC (I) 


INC ( J) 


CLDEN(I) 


Altitude of bottom surface of 
cloud (in km.) 

Edge lengths of I*'* 1 cloud (in km.) 

Angle between the x-axis of the 
cloud and that of the general 
co-ordinate system (measured 
clockwise in degrees looking up 
from the ground) . 

Indicator of the atmospheric con- 
stituent appropriate to I L cloud 
[e-. g. , INC (1) = 3 states that cloud 
number 1 is composed of the species 
specified by INDIC(3)]. 

Density of cloud-specific scatters 
of the type labeled by INC(I) (in cm ^) . 
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IRCLD , I REC , RECANG , HR , CZO , RPHI , RX , RY , SZO , I AP AT , NCON 
(212 , 7E10 .4,212) 

IRCLD = Receiver extent indicator 

= 0 Receiver with infinitesimal area 
= J>0 Receiver extends over an entire face 
of cloud J. 

IREC Specifies which face of the cloud is to be 

the receiver. 

= 1 + X Face of cloud IRCLD 

= 2 - X Face of cloud IRCLD 

= 3 + Y Face of the cloud IRCLD 

= 4 - Y Face of cloud IRCLD 

= 5 + Z Face of cloud IRCLD 

= 6 - Z Face of Cloud IRCLD 

< 0 Receiver points out of cloud IRCLD. 
j IREC | specifies which face as above. 

RECANG = Apex angle of cone of receiver field 
of view in milliradians . 

HR = Height of receiver (ignored if receiver is 

a cloud face) 

CZO = Zenith cosine of receiver pointing direction 

(ignored if receiver is a cloud face) 

RPHI = Azimuthal angle of plane containing receiver 

field of view and the vertical (ignored if 
receiver is a cloud face) 

RX,RY = X,Y co-ordinates of receiver (ignored if 
receiver is a cloud face) 

SZO = Sine of receiver zenith angle 

= 0 . or blank - SZO is computed from CZO 
f 0. CZO,, computed from SZO (ignored if 
receiver is a cloud face) 

IAPAT = Index of albedo pattern for which the intensity 
is to be mapped within the receiver field of 
view . 

DEFAULT VALUE- IAPAT=1 (Max=50) 

NCON = Number of THETA BINS=Number of PHI bins 

for mapping of response within the receiver 
field of view. 

DEFAULT=20 (Max=20) 
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IRPAT , IRBGD , XBND , FN , WSPD , WCO , WC1 , WS IG2 
(212 , 6E11 . 4) 

= Specification of reflectance type map 
= 0 Uniform plane (= Background) 

=+l Plane differs from background for X>XBND 

=-l Plane differs from background for X<XBND 

= 2 Inhomogeneties are rectangles set upon 
uniform background plane. 

= Specification of character of background 
= 0 Background is Lambertian 

= 1 Background is rough Fresnel 

= Boundary line for half-plane pattern 
= X-co-ordinate of line parallel to Y-axis 
= Index of refraction for Fresnel reflectance 
DEFAULT= 1.338 

= Wind speed (m/sec.) for rough Fresnel 
= Parameters used to construct variance, WSIG2 
If both = 0 (or Blank), use default values: 

WC0= .0015 
WC1= 2 . 54E- 03 

= Variance for Gaussian slope distribution for 
rough Fresnel reflectance 

Non- zero or non-blank entry takes precedence 
over value computed from WC0,WC1, and WSPD. 

If IRPAT=2, the following cards are required. 

NSPOT 

( 12 ) 

NSPOT = Number of non- overlapping rectangular inhomo- 
geneities (maximum NSPOT=20) 

CNTRX(l) , CNTRY (1) ,EDGEX(1) ,EDGEY (1) 

(4E11.4) 

CNTRX (NSPOT) , CNTRY (NSPOT) , . . . 

CNTRX (J) , CNTRY (J) = X-Y co-ordinates of rectangle J 
EDGEX (J) ,EDGEY (J) = X,Y edge lengths of rectangle J 


IRPAT 

IRBGD 

XBND 

FN 

WSPD 

WC0,WC1 

WSIG2 
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MIE TAPE FORMATS 


MONODISPERSE 

/ REFACT, ALPHA, QSCA , COSBAR , QEXT 

I = 1 ’ n 

* J=1 ,NK (K) { XA , AA1 , AA2 ,S33(J,K,1) ,S34(J,K,1) ,S11(J,K) ,S12(J,K,1) 

where 


n is computed from A1(K) and STEP(K) read in on cards: 

n = INT (A1 (K) (STEP (K) + . 5) 

REFACT = complex index of refraction 

ALPHA = size parameter a = ud/X, d = diameter, 

X = wavelength 

QSCA = scattering cross section normalized by geo- 
metric cross section 

QEXT = extinction cross section normalized by geo- 
metric cross section 

2 d^ - 8 

eg a SCA ( - Cm ^ = ^SCA X 17 "4" X 10 

with d in micrometers 

COSBAR - not used 

NK(K) = number of scattering angles for species k 

XA = scattering angle (degrees) 


c * 


94a 


533 

511 

512 

534 = 


Scattering matrix elements for Jth scattering 
angle kth species 

-S34 as used by this code 


POLYDISPERSE 

NT, ST [Read only if ITAPE^O; for ITAPE=0, values are 
taken from card input.] 

( REFACT , ALPHA , QSCA , COSBAR , QEXT 
1=1, NT) 

(j=l,NK(K) { XA , AA, BB , XX3 , XX4 , YY , XX2 
NT = number of steps in size 

ST = dimensionless step size = ^ x id where Sd is the 
step size in diameter (microns) 

REFACT, ALPHA, QSCA, COSBAR, QEXT as described under 
MONODISPERSE input 


XA = scattering angle (degrees) 


AA,BB 
XX3 = 
XX4 = 
YY = 
XX 2 = 


- not used 



Scattering matrix elements for size 
parameter ALPHA, scattering angle XA 
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SECTION 6.0 

CTRANS PROGRAM LISTING 
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